Numerical recovering a density by BC-method 



In this paper we develop numerical algorithm for solving inverse problem for the wave 
equation using Boundary Control method. The results of numerical experiments are represented. 
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1 Introduction 

The Boundary Control method is one of the most natural method for solving 
multidimensional inverse dynamical problems. The method was proposed by Belishev in 1986 (see 
[B], [Bl] and publications cited there). Numerical testing results are represented in [BG]. 

The principle question in the BC-method is the controllability. For scalar equations like 
wave equation the approximate controllability is valid. It means any function (p in bounded 

domain Q (say from may be arbitrary close approximated by a wave u^{.,T), induced by 

some boundary source / (control). For arbitrary <p the equation u^{.,T) = q> w.r.t. / cannot be 
solved when coefficients of the (wave) equation are unknown. The remarkable fact is that for 
harmonic g) one can find a control which provides the equality u\.,T) = (p with arbitrary 
accuracy using only data of the inverse problem. But in practice we have only finite number of 
controls. Can one provide proper accuracy at that and how many controls is takes? In this paper we 
try to answer these questions numerically when reconstructing a density in the unit disk in . We 
take T more than optical radius of Q. It makes our problem easier. 

We prove also approximate i/' -controllability. We use this to develop a numerical 
algorithm for solving inverse problem for the wave equation with unknown density. We also 
demonstrate the results of numerical experiments. 

Let Q be a bounded one-connected domain in R" (n = 2,3) with a smooth boundary F. 
Denote by the solution (wave) to the initial boundary problem for the wave equation 



where p{x)>0 is a smooth density, m„ is the normal derivative (control). Here Z^(rx0,2r]) is 
real (as all Hilbert spaces in this paper). The map R^^ : -F^ lI^T x 0,2r]) defined by 



yC»M„ - Am = /« Q X (0,2 T) 



(1) 
(2) 
(3) 




\t>T 



0}, 




is classical when /gM, 



Consider inverse problem: to reconstruct p from R under assumption 



T >T* = sup^^^ dist [x, r), (4) 
where the distance is tal<ing with respect to the metric yjp{x)\dx\ . This assumption provides that 
the waves induced by various controls fill up the closure of Q at the final moment T . Note, that if 
0<r<r*then data {f,R^^f),f gF^ of the the inverse problem do not include information 

about the density on the set Q\Q^, where 

= {jc e Q I < r} 

is the set filled by waves up to the moment T . In our case =Q. In what follows the fixed final 
moment T>T* will not be mentioned at all notations. 

Our approach to this problem is based on the BC-method and close to approach of [P]. 

2 Bilinear forms 

Introduce two symmetric bilinear forms 

Vf,g] = \^^{.,Ty{.,T)dx, (5) 

The following relations between these forms and response operator are the base of BC-method. 
We obtain these relations for the convenience of riders (for more details see [B]). We use the 
notations 

(., = t)±u{.,2T- t))l2, 
ilf){.,t) = \/{.,s)ds, ?e0,2r], 
dL is the volume form of Fx 0,r]. 

Proposition 1 For any controls /, g e M the equalities 

= l^^^l{Rg)Jf-gJRn<^, (6) 

are valid. 

Proof. For any solution v to the wave equation the equality 

p(vuf - M^v,), = div(vVu^ - u^Vv) (8) 

holds. Substituting uf for v and integrating over cylinder QxO,r], we get 

Ipu'i., T)u{{., T)dx = l^^^lulf - u^{ul)„W. 

Taking into account that u' =ul and u ^ = lu ' , we get (6). 

Consider (potential) bilinear form [f,g'\p. For any solution v to the wave equation the 
equality 

p{v,u{),+{Wv,Wu'),=div{v;^u^ + u{Wv). 
holds. Substituting t/f for v and integrating over [Q x 0,r], we get 



Remark 2 Analogously to (7) one can get representation of the kinetic form 

def r f) r) 

[/, g], = \^{x)u! {x, T)uf {x, T)dx = ^^^[/ - {Rg)_ +g_- RfW. 0) 



Remark 3 ft may be shown (6),(7),(9) do not depend on values of controls on the set 
rxT,2T]. 



3 Boundary control and density reconstruction 

Consider the boundary control problem 

u^{.,T) = q?^L'{n), / = ? (10) 

For any T>0 the linear variety H = {u^(., T)\f e M } is dense in l}(fl) [B]. For sufficiently large 

T>T* and under some geometrical assumption the equation is solvable in F (not uniquely) 
[BLR]. 

We consider the case when <pgH\Q.), where H\Q.) is the real Hilbert space with the 

norm 

= (^jou^dx + 1^ I Vu f dxf'^. 
Show that the set H is dense in H\Q). 

Theorem 4 The orthogonal complement to H in H\D.) is {0}. 
Proof. Let v be the solution to the initial boundary value problem 

yCn;„-Av = 0,mQx(0,r), 

lrxo,r]~ 

where cp e . Since (p e //'(Q) we have 

V e C([0, T\H\m, V, e C([0, r];/f'(Q)), v„ e C([0, r];Z^(Q)) 
[LM]. Denote by u the wave m = e M . By standard way 

J^("v, - vuX; T)dx = j^^g^j(wv„ - vuJdL =^ 

T)(pdx - £Am^(., T)(pdx = l^^^^v{f„ - f)dZ ^ 
jjM^ i., T)(pdx + \jyu' (., n V(p)dx = Iqfi., T)dL + v(/, - f)dL 

Thus we have 



{uf,(p) 1 = f {vJ, + vf)dZ. 
^^hUo) JrxOTi 



v{.,t) 



'(n) Jrxo.T] 

Therefore implies v |p^oyj= 0. The even extension of v w.r.t. t = T 

v{.,t), t&0,T] 
-v{.,2r-t), tGT,2T] 
is the function from C([0, T];H^{Ciy) which satisfies 

pv„ - Av = 0, in (0,2 T), 

^ lrxo,2r]~ 0' ^« lrxo,2r]~ 0- 

As in [B] using Tataru's theorem [T] this implies v = in the domain bounded by characteristics 

t = dist(x,r) and t = 2T -dist(x, T) .Thus (p = 0. 

Let (p be an arbitrary smooth harmonic function in QuF. Consider the functional 

^■.H\rxO,T])^R: 

def _ 

^{f) = \^\Vu\.,T)-V(pfdx 

= f \Vu^ f i.,T)dx-2\ (yu^{.,T)y(p)dx+ \ \V(pfdx 

JQ JQ Jfi 

= [/,/], - 2l{Rf){x,T)(p^{x)dr + Icpcp^dT. 

Note, that O is completely determined by response operator. The following propositions is the 
base of our approach to reconstruct p . 

Proposition 5 For any smooth harmonic function cp and s>0 there is a control / e M 
such that 



Inequality (11) implies 



0(f)+\\Rf(.,T)-^\f^,^^^<s. (11) 



\x,T)-(p\f, <C£, 



with some constant C, does not depending on / and ^. 

Proof. The first statement follows from the density H in H\D.) and the definition of O. 
The second one follows from Friederischs' inequality. 

Thus one can control the closeness u^(., T) to cp in //' from the boundary. 
The reconstruction of the density may be fulfilled by the following scheme. 

1. For any smooth harmonic functions q) one can find the control e M such that 

^\{f,p)+\\Rf<p{-J)-9\\ji^^^'^^ and therefore \\u''{x,T)-(p\\^^^^^<Cs with arbitrarily small £>0 

2. Substituting /^^ for and /^^ for in (5) we get approximate equality 

[/'(^ Vi V2 {x)dx « [/^j , /pj , (12) 
where (p^,(p2 are arbitrary smooth harmonic functions. Since the linear span of all products (p^cpj is 
dense in l}{fl) then one can use (12) to find p. When numerical solving the inverse problem we 



use also the a priori limitations for p . Emphasize that both procedures are linear. 



4 Discrete inverse problem 

Project the forward problem (1-3) onto a finite dimensional space that spans standard 
continuous piecewise basic functions i//^{x),n = l,..., N of the Finite Element Method (A^ is the 
number of nodes). Then (1-3) reduces to the Cauchy problem for a linear system of ordinary 
differential equations. 

The projection u{/ of the wave is expended into the finite sum 

N 
n=l 

where x„ are nodes. Vector-function U-^ is the solution of the Cauchy problem for ordinary 
differential equations with constant coefficients : 

MU„+KU = Gf, G({t)=\xi^j{.,t)dr (13) 

u{o)=uXo)=o. 

where 

are mass matrix and stiffness matrix accordingly. 

As in the differential case the following representations 

def .J 

U,gf ={U',MWtT) = I [iIG',U'J- {GlJU'mdt, (14) 

JO 

Vf,gt =iU' ,KU^){T) = lV{Gf ,Ur) + {Gl,U{m)dt. (15) 
holds. Emphasize that right hand sides of (14), (15) depend on l.^.o^n only, where r, is the set 
of boundary nodes. 

Numerical experiments was made for unit disk Q with continuous piecewise constant 
model of the density and the controls 

fja{^,t)=r{t-jM)qSx\ 7 = 1,..., N, = TIM, a = \,..., N,. (16) 

Here r{t) is Ricker's impulse (fig.l), are continuous piecewise "linear" functions on the 

boundary, qaixp) = 5^p,\xp \=\,N^ is the number of boundary nodes. In what follows the basic 

controls (16) are numbered by one index i = 1,..., iV^ xiV^ . 



1. Ricker wavelet 



The matrixes 

C,=Vf„fjf, P,=Vf,Jjt, \<i,j<N,^N, 
and values U^'(j)\^^ were used as the data of inverse problem. The reconstruction of p was 

fulfilled using a scheme close to described above for differential inverse problem. 

1. Harmonic mesh functions. Define harmonic mesh functions (p^,a = l,...N^ as solutions 
of equations 

K(Pa=K, a = l,...,N,. 
where linearly independent vectors satisfy conditions of solvability 

Xi^„(x,) = 0,a = l,...^,, 

The final harmonic mesh function is : 

VN^^+ii^j) = 1, 7 = I,--- A, K(p^^^, = 0. 
The q>^ is the mesh analogue of a harmonic function satisfying to the Neumann condition. 

2. The control problem. Formally substituting (p^ for U* , and / for g in (15) we get 
equations 

U. Jt= (U'- , Z„ ){T), i = l,...JM„xN, (17) 
w.r.t. /. To find a control, which gives (t) » (p^ we have to use also equalities 

[/^(x,,r) = ^„(x,), x,eT,. (18) 
Denote by the normal solution to the system (17),(18). The control provided a well accuracy 
of equality i7^(r)«^^ in each our numerical experiment (relative error in P is about 10 **). 

3. Reconstruction. Substituting / for f^,g)^ for U' (T), g for f^, and g?^ for in (14) 

we get 

((P„,M(P,) = ILJ.r, a,/? = 1,... 
Let p,^ be the value of p in the k"' triangle, k = l,...JC. Then we get the linear algebraic system 
w.r.t. p^. : 

tj^,(P>il ¥i{xVj{x)dx = UaJpf, a,p = l,..., N, + l. (19) 
This system was sometimes ill-conditioned. Therefore we used natural a priori limitations for 



values pi^ and optimization algorithms to reconstruction. Below the results of numerical 
experiments are represented, S means relative error in f. 




4. Sample 2. Free inner boundary 5. Reconstruction of sample 2, 5 — \Wo 




6. Sample 3. Waveguides 



7. Reconstruction of sample 3, 5 = 2, 5% 



8. Sample 4. Folds 



9. Reconstruction of sample 4, (5 = 3% 
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